Lower-Bounded Analysis

0.1 Setup

library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)
library(ggthemes)
library(viridis)
library(brew)
library(latex2exp)
library(scico)
library(DescTools)
library(brms)
library(ggh4x)

# Change this to the directory you want the figures to be saved to
figure_path = "./figures/"

# Change this to the directory that contains the full_sim_result.RDS file
data_path = "./"

# This is used for brms model fitting
NCORES = 12

sum_coding <- function(x, lvls = levels(x)) {
  # codes the first category with -1
  nlvls <- length(lvls)
  stopifnot(nlvls > 1)
  cont <- diag(nlvls)[, -nlvls, drop = FALSE]
  cont[nlvls, ] <- -1
  cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
  colnames(cont) <- lvls[-1]
  x <- factor(x, levels = lvls)
  contrasts(x) <- cont
  x
}

0.2 Preprocess

result_df <- readRDS(paste0(data_path, "full_sim_result.RDS"))
result_df <- result_df %>%
filter(fit_family != "inverse.gaussian",
       formula != "y ~ x + z2" & formula != "y ~ x + z1 + z2 + z4",
        divergents <= 10,
        rhat <= 1.01,
        ess_bulk > 400,
        ess_tail > 400) %>%
  # pretty format all the labels for figure printing  
mutate(
  fit_link_transformed = case_when(
    fit_family == "lognormal" ~ "Log",
    fit_family == "softplusnormal" ~ "Softplus",
    fit_link == "log" ~ "Log",
    fit_link == "softplus" ~ "Softplus",
    fit_link == "identity" ~ "Identity",
    TRUE ~ fit_link
  )
) %>%
mutate(
  data_link_transformed = case_when(
    data_family == "lognormal" ~ "Log",
    data_family == "softplusnormal" ~ "Softplus",
    data_link == "log" ~ "Log",
    data_link == "softplus" ~ "Softplus",
    TRUE ~ data_link
  )
) %>%
mutate(
  fit_family_transformed = case_when(
  fit_family == "lognormal" |
  fit_family == "softplusnormal" ~  "trans. Normal",
  fit_family == "gaussian" ~ "Normal",
  fit_family == "gompertz" ~ "Gompertz",
  fit_family == "gamma" ~ "Gamma",
  fit_family == "weibull" ~ "Weibull",
  fit_family == "frechet" ~ "Fréchet",
  fit_family == "betaprime" ~ "Beta prime",
  TRUE ~ fit_family
  )
) %>%
mutate(
  data_family_transformed = case_when(
  data_family == "lognormal" |
  data_family == "softplusnormal" ~  "trans. Normal",
  data_family == "gaussian" ~ "Normal",
  data_family == "gompertz" ~ "Gompertz",
  data_family == "gamma" ~ "Gamma",
  data_family == "weibull" ~ "Weibull",
  data_family == "frechet" ~ "Fréchet",
  data_family == "betaprime" ~ "Beta prime",
  TRUE ~ data_family
  )
) %>%
  mutate(shape = case_when(
    shape == "symmetric" ~ "Thin Tail",
    shape == "asymmetric" ~ "Heavy Tail",
    shape == "ramp" ~ "Ramp"
  )) %>%
  # recreate the dataset id so it is truly unique
  mutate(dataset_id = as.numeric(as.factor(paste0(
    dataset_seed,
    data_family_transformed,
    data_link_transformed,
    x_y_coef,
    y_intercept,
    sigma_y
  )))) %>%
  mutate(rmse_loo = abs(rmse_loo),
         rmse_newdata = abs(rmse_newdata),
         abs_bias = abs(bias)) %>%
  mutate(formula = factor(formula, levels = c("y ~ x + z1 + z2",
                                                 "y ~ x + z1",
                                                 "y ~ x + z1 + z2 + z3")),
         shape = factor(shape, levels  = c("Thin Tail",
                                           "Heavy Tail",
                                           "Ramp")),
         data_family_transformed = factor(data_family_transformed, levels =
                                            c("Gamma",
                                              "Weibull",
                                              "Fréchet",
                                              "Beta prime",
                                              "Gompertz",
                                              "trans. Normal")),
         fit_family_transformed = factor(fit_family_transformed, levels =
                                            c("Gamma",
                                              "Weibull",
                                              "Fréchet",
                                              "Beta prime",
                                              "Gompertz",
                                              "trans. Normal",
                                              "Normal")))

0.3 RMSE & Bias

cols <- c("dataset_id", "formula", "shape", "fit_link_transformed")

df <- result_df %>%
  filter(fit_link_transformed == data_link_transformed) %>%
  mutate(
     fit_link_transformed = factor(fit_link_transformed, levels = 
                                         c("Log",
                                           "Softplus")),
   data_link_transformed = factor(data_link_transformed, levels = 
                                          c("Log",
                                           "Softplus"))
   )

# full RMSE results
p_rmse <- df %>%
  ggplot(aes(y = fit_family_transformed, x = rmse_s)) +
  coord_cartesian(xlim = c(0, 1)) +
  scale_x_continuous(labels = c("0", "0.5", "1"), breaks = c(0, 0.5, 1)) +
  geom_boxplot() +
  stat_summary(fun="mean", shape = 4, size = 0.3) +
  facet_grid(data_link_transformed + shape ~ data_family_transformed) +
  xlab(TeX("$RMSE(beta_{xy})$")) +
  ylab("Fit likelihood") +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        axis.ticks.y = element_blank())
  
p_rmse

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

ggsave(paste0(figure_path, "positive_rmse.pdf"), width = 210, height = (297/4) * 3, units = "mm", useDingbats = TRUE)

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

# RMSE log showcase
scales <- list(
  scale_x_continuous(limits = c(0, 0.2), breaks = c(0, 0.1, 0.2), labels = c("0", "0.1", "0.2")),
  scale_x_continuous(limits = c(0, 0.3), breaks = c(0, 0.1, 0.2, 0.3), labels = c("0", "0.1", "0.2", "0.3")),
  scale_x_continuous(limits = c(0, 0.5), breaks = c(0, 0.2, 0.4), labels = c("0", "0.2", "0.4"))
)
p_rmse_pointrange <- df %>%
  filter(data_link_transformed == "Log") %>%
  group_by(data_family_transformed, fit_family_transformed, shape) %>% 
  summarise(median_rmse_s = median(rmse_s),
            min_rmse = quantile(rmse_s, probs = c(0.025)),
            max_rmse = quantile(rmse_s, probs = c(0.975))) %>%
  ggplot(
    aes(y = data_family_transformed,
        x = median_rmse_s,
        xmin = min_rmse,
        xmax = max_rmse,
        color = fit_family_transformed)) +
  geom_pointrange(size = 0.18, position = position_dodge(width = 0.8)) +
  scale_color_scico_d(palette = "batlow", end = 0.8) +
  facet_grid(~ shape, scales = "free_x") +
  facetted_pos_scales(x = scales) +
  xlab(TeX("$RMSE(beta_{xy})$")) +
  ylab("Data likelihood") +
  labs(color = "Fit likelihood") +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        axis.ticks.y = element_blank())

## `summarise()` has grouped output by 'data_family_transformed',
## 'fit_family_transformed'. You can override using the `.groups` argument.

p_rmse_pointrange

ggsave(paste0(figure_path,"positive_rmse_showcase_pointrange.pdf"), width = 210, height = (297/4)*1.4, units = "mm", useDingbats = TRUE)

# RMSE softplus showcase
p_rmse_softplus_showcase <- df %>%
  filter(data_link_transformed == "Softplus") %>%
  ggplot(aes(y = fit_family_transformed, x = rmse_s)) +
  coord_cartesian(xlim = c(0, 0.75)) +
  scale_x_continuous(labels = c("0", "0.75"), breaks = c(0, 0.75)) +
  geom_boxplot() +
  stat_summary(fun="mean", shape = 4, size = 0.3) +
  facet_grid(shape ~ data_family_transformed) +
  xlab(TeX("$RMSE(beta_{xy})$")) +
  ylab("Fit likelihood") +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        axis.ticks.y = element_blank())

p_rmse_softplus_showcase

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

ggsave(paste0(figure_path,"positive_rmse_softplus_showcase.pdf"), width = 210, height = (297/4)*1.5, units = "mm", useDingbats = TRUE)

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

# full absolute bias results
p_bias <- df %>%
  ggplot(aes(y = fit_family_transformed, x = abs_bias)) +
  coord_cartesian(xlim = c(0, 0.5)) +
  scale_x_continuous(labels = c("0", "0.25", "0.5"), breaks = c(0, 0.25, 0.5)) +
  geom_boxplot() +
  stat_summary(fun="mean", shape = 4, size = 0.3) +
  facet_grid(data_link_transformed + shape ~ data_family_transformed) +
  xlab(TeX("$abs(bias(beta_{xy}))$")) +
  ylab("Fit likelihood") +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        axis.ticks.y = element_blank())
  
p_bias

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

ggsave(paste0(figure_path, "positive_bias.pdf"), width = 210, height = (297/4) * 3, units = "mm", useDingbats = TRUE)

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

# absolute bias showcase
p_bias_showcase <- df %>%
  filter(data_link_transformed == "Log") %>%
  ggplot(aes(y = fit_family_transformed, x = abs_bias)) +
  coord_cartesian(xlim = c(0, 0.25)) +
  scale_x_continuous(labels = c("0", "0.25"), breaks = c(0, 0.25)) +
  geom_boxplot() +
  stat_summary(fun="mean", shape = 4, size = 0.3) +
  facet_grid(shape ~ data_family_transformed) +
  xlab(TeX("$abs(bias)(beta_{xy})$")) +
  ylab("Fit likelihood") +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        axis.ticks.y = element_blank())

p_bias_showcase

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

ggsave(paste0(figure_path,"positive_bias_showcase.pdf"), width = 210, height = (297/4)*1.5, units = "mm", useDingbats = TRUE)

## Warning: Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).
## Removed 7 rows containing missing values (`geom_segment()`).

0.4 ROC

calibration_df <- result_df %>%
  # Calculate significance of the .
  mutate(x_y_coef_cat = ifelse(x_y_coef == 0, 0, 1),
         sig95 = pq_0.025 > 0 | pq_0.975 < 0,
         sig90 = pq_0.05 > 0 | pq_0.95 < 0,
         sig80 = pq_0.1 > 0 | pq_0.9 < 0,
         sig50 = pq_0.25 > 0 | pq_0.75 < 0)

# We add the normal identity model to all fit links for easier visual comparison
calibration_df <- calibration_df %>%
   mutate(
    fit_family_transformed = case_when(
    (fit_family_transformed == "Normal" & fit_link_transformed == "Identity") ~ "Normal Identity",
    TRUE ~ fit_family_transformed)) %>%
  mutate(
    fit_link_transformed = case_when(
    (fit_family_transformed == "Normal Identity") ~ "Log",
    TRUE ~ fit_link_transformed))
data1 <- calibration_df %>%
  filter(fit_family_transformed == "Normal Identity") %>%
  mutate(fit_link_transformed = "Softplus")

calibration_df <- rbind(calibration_df, data1)

calibration_df = calibration_df %>%
  mutate(fit_link_transformed = factor(fit_link_transformed, levels = 
                                         c("Log",
                                           "Softplus")))

cols = c("formula", "fit_link_transformed", "fit_family_transformed", "shape", "data_link_transformed", "data_family_transformed")

TPR <- calibration_df %>%
filter(x_y_coef != 0) %>%
group_by(across(all_of(cols))) %>%
summarise(TPR_95 = mean(sig95),
          TPR_90 = mean(sig90),
          TPR_80 = mean(sig80),
          TPR_50 = mean(sig50),
          TPR_0 = 0,
          TPR_1 = 1)

## `summarise()` has grouped output by 'formula', 'fit_link_transformed',
## 'fit_family_transformed', 'shape', 'data_link_transformed'. You can override
## using the `.groups` argument.

FPR <- calibration_df %>%
filter(x_y_coef == 0) %>%
group_by(across(all_of(cols))) %>%
summarise(FPR_95 = mean(sig95),
          FPR_90 = mean(sig90),
          FPR_80 = mean(sig80),
          FPR_50 = mean(sig50),
          FPR_0 = 0,
          FPR_1 = 1)

## `summarise()` has grouped output by 'formula', 'fit_link_transformed',
## 'fit_family_transformed', 'shape', 'data_link_transformed'. You can override
## using the `.groups` argument.

df = merge(TPR, FPR)
auc <- vector(mode = "numeric", length = nrow(df))
for (i in seq_len(nrow(df))) {
  auc[[i]] <- AUC(
    x = c(df$FPR_0[[i]], df$FPR_95[[i]], df$FPR_90[[i]], df$FPR_80[[i]], df$FPR_50[[i]],df$FPR_1[[i]]),
    y = c(df$TPR_0[[i]], df$TPR_95[[i]], df$TPR_90[[i]], df$TPR_80[[i]], df$TPR_50[[i]], df$TPR_1[[i]])
  )
}

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

df$auc <- auc
df <- df %>%
  mutate(formula = sum_coding(
    factor(formula,
           levels = c("y ~ x + z1 + z2",
                      "y ~ x + z1",
                      "y ~ x + z1 + z2 + z3")))) %>%
  mutate(data_link_transformed = sum_coding(
    factor(data_link_transformed, levels = c("Log",  "Softplus"))
  )) %>%
  mutate(fit_family_transformed = sum_coding(
    factor(fit_family_transformed,
           levels = c("Gamma",
                      "Weibull",
                      "Beta prime",
                      "trans. Normal",
                      "Fréchet",
                      "Gompertz",
                      "Normal",
                      "Normal Identity")))) %>%
  mutate(shape = sum_coding(
    factor(shape, levels = c("Thin Tail", "Heavy Tail", "Ramp")))
  )

m1 = brm(auc ~ 1 +
         fit_link_transformed + fit_family_transformed +
         data_link_transformed + data_family_transformed +
         shape + formula +
         fit_link_transformed * data_link_transformed +
         fit_family_transformed * data_family_transformed,
         data = df, family = Beta(), cores = 4, warmup = 1000, iter = 3500,file = paste0(data_path, "m_auc.RDS"))
summary(m1)

##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: auc ~ 1 + fit_link_transformed + fit_family_transformed + data_link_transformed + data_family_transformed + shape + formula + fit_link_transformed * data_link_transformed + fit_family_transformed * data_family_transformed 
##    Data: df (Number of observations: 1725) 
##   Draws: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup draws = 10000
## 
## Population-Level Effects: 
##                                                                          Estimate
## Intercept                                                                    2.02
## fit_link_transformedSoftplus                                                -0.18
## fit_family_transformedWeibull                                                0.07
## fit_family_transformedBetaprime                                              0.04
## fit_family_transformedtrans.Normal                                           0.32
## fit_family_transformedFréchet                                               -0.74
## fit_family_transformedGompertz                                              -0.23
## fit_family_transformedNormal                                                 0.01
## fit_family_transformedNormalIdentity                                         0.15
## data_link_transformedSoftplus                                               -0.01
## data_family_transformedWeibull                                               0.14
## data_family_transformedFréchet                                               0.25
## data_family_transformedBetaprime                                             0.23
## data_family_transformedGompertz                                             -0.09
## data_family_transformedtrans.Normal                                         -0.32
## shapeHeavyTail                                                               0.02
## shapeRamp                                                                    0.06
## formulay~xPz1                                                                0.03
## formulay~xPz1Pz2Pz3                                                         -0.24
## fit_link_transformedSoftplus:data_link_transformedSoftplus                   0.20
## fit_family_transformedWeibull:data_family_transformedWeibull                 0.27
## fit_family_transformedBetaprime:data_family_transformedWeibull              -0.26
## fit_family_transformedtrans.Normal:data_family_transformedWeibull            0.00
## fit_family_transformedFréchet:data_family_transformedWeibull                -0.45
## fit_family_transformedGompertz:data_family_transformedWeibull                0.25
## fit_family_transformedNormal:data_family_transformedWeibull                  0.04
## fit_family_transformedNormalIdentity:data_family_transformedWeibull          0.10
## fit_family_transformedWeibull:data_family_transformedFréchet                -0.85
## fit_family_transformedBetaprime:data_family_transformedFréchet               0.75
## fit_family_transformedtrans.Normal:data_family_transformedFréchet            0.10
## fit_family_transformedFréchet:data_family_transformedFréchet                 1.50
## fit_family_transformedGompertz:data_family_transformedFréchet               -0.66
## fit_family_transformedNormal:data_family_transformedFréchet                 -0.17
## fit_family_transformedNormalIdentity:data_family_transformedFréchet         -0.24
## fit_family_transformedWeibull:data_family_transformedBetaprime              -0.34
## fit_family_transformedBetaprime:data_family_transformedBetaprime             0.45
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime          0.19
## fit_family_transformedFréchet:data_family_transformedBetaprime               0.46
## fit_family_transformedGompertz:data_family_transformedBetaprime             -0.37
## fit_family_transformedNormal:data_family_transformedBetaprime               -0.03
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime       -0.08
## fit_family_transformedWeibull:data_family_transformedGompertz                0.48
## fit_family_transformedBetaprime:data_family_transformedGompertz             -0.83
## fit_family_transformedtrans.Normal:data_family_transformedGompertz          -0.11
## fit_family_transformedFréchet:data_family_transformedGompertz               -0.75
## fit_family_transformedGompertz:data_family_transformedGompertz               0.71
## fit_family_transformedNormal:data_family_transformedGompertz                 0.30
## fit_family_transformedNormalIdentity:data_family_transformedGompertz         0.25
## fit_family_transformedWeibull:data_family_transformedtrans.Normal           -0.18
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal          0.17
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal       0.01
## fit_family_transformedFréchet:data_family_transformedtrans.Normal            0.37
## fit_family_transformedGompertz:data_family_transformedtrans.Normal          -0.06
## fit_family_transformedNormal:data_family_transformedtrans.Normal            -0.00
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal    -0.06
##                                                                          Est.Error
## Intercept                                                                     0.03
## fit_link_transformedSoftplus                                                  0.03
## fit_family_transformedWeibull                                                 0.08
## fit_family_transformedBetaprime                                               0.08
## fit_family_transformedtrans.Normal                                            0.09
## fit_family_transformedFréchet                                                 0.07
## fit_family_transformedGompertz                                                0.08
## fit_family_transformedNormal                                                  0.08
## fit_family_transformedNormalIdentity                                          0.08
## data_link_transformedSoftplus                                                 0.02
## data_family_transformedWeibull                                                0.05
## data_family_transformedFréchet                                                0.05
## data_family_transformedBetaprime                                              0.05
## data_family_transformedGompertz                                               0.04
## data_family_transformedtrans.Normal                                           0.04
## shapeHeavyTail                                                                0.02
## shapeRamp                                                                     0.02
## formulay~xPz1                                                                 0.02
## formulay~xPz1Pz2Pz3                                                           0.02
## fit_link_transformedSoftplus:data_link_transformedSoftplus                    0.03
## fit_family_transformedWeibull:data_family_transformedWeibull                  0.13
## fit_family_transformedBetaprime:data_family_transformedWeibull                0.11
## fit_family_transformedtrans.Normal:data_family_transformedWeibull             0.13
## fit_family_transformedFréchet:data_family_transformedWeibull                  0.10
## fit_family_transformedGompertz:data_family_transformedWeibull                 0.11
## fit_family_transformedNormal:data_family_transformedWeibull                   0.12
## fit_family_transformedNormalIdentity:data_family_transformedWeibull           0.12
## fit_family_transformedWeibull:data_family_transformedFréchet                  0.11
## fit_family_transformedBetaprime:data_family_transformedFréchet                0.14
## fit_family_transformedtrans.Normal:data_family_transformedFréchet             0.14
## fit_family_transformedFréchet:data_family_transformedFréchet                  0.13
## fit_family_transformedGompertz:data_family_transformedFréchet                 0.11
## fit_family_transformedNormal:data_family_transformedFréchet                   0.12
## fit_family_transformedNormalIdentity:data_family_transformedFréchet           0.12
## fit_family_transformedWeibull:data_family_transformedBetaprime                0.12
## fit_family_transformedBetaprime:data_family_transformedBetaprime              0.13
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime           0.14
## fit_family_transformedFréchet:data_family_transformedBetaprime                0.11
## fit_family_transformedGompertz:data_family_transformedBetaprime               0.11
## fit_family_transformedNormal:data_family_transformedBetaprime                 0.12
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime         0.12
## fit_family_transformedWeibull:data_family_transformedGompertz                 0.12
## fit_family_transformedBetaprime:data_family_transformedGompertz               0.11
## fit_family_transformedtrans.Normal:data_family_transformedGompertz            0.12
## fit_family_transformedFréchet:data_family_transformedGompertz                 0.09
## fit_family_transformedGompertz:data_family_transformedGompertz                0.12
## fit_family_transformedNormal:data_family_transformedGompertz                  0.12
## fit_family_transformedNormalIdentity:data_family_transformedGompertz          0.12
## fit_family_transformedWeibull:data_family_transformedtrans.Normal             0.11
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal           0.11
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal        0.12
## fit_family_transformedFréchet:data_family_transformedtrans.Normal             0.10
## fit_family_transformedGompertz:data_family_transformedtrans.Normal            0.10
## fit_family_transformedNormal:data_family_transformedtrans.Normal              0.11
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal      0.11
##                                                                          l-95% CI
## Intercept                                                                    1.95
## fit_link_transformedSoftplus                                                -0.23
## fit_family_transformedWeibull                                               -0.09
## fit_family_transformedBetaprime                                             -0.12
## fit_family_transformedtrans.Normal                                           0.15
## fit_family_transformedFréchet                                               -0.88
## fit_family_transformedGompertz                                              -0.38
## fit_family_transformedNormal                                                -0.15
## fit_family_transformedNormalIdentity                                        -0.01
## data_link_transformedSoftplus                                               -0.04
## data_family_transformedWeibull                                               0.05
## data_family_transformedFréchet                                               0.16
## data_family_transformedBetaprime                                             0.14
## data_family_transformedGompertz                                             -0.18
## data_family_transformedtrans.Normal                                         -0.41
## shapeHeavyTail                                                              -0.02
## shapeRamp                                                                    0.03
## formulay~xPz1                                                               -0.01
## formulay~xPz1Pz2Pz3                                                         -0.27
## fit_link_transformedSoftplus:data_link_transformedSoftplus                   0.15
## fit_family_transformedWeibull:data_family_transformedWeibull                 0.03
## fit_family_transformedBetaprime:data_family_transformedWeibull              -0.48
## fit_family_transformedtrans.Normal:data_family_transformedWeibull           -0.25
## fit_family_transformedFréchet:data_family_transformedWeibull                -0.64
## fit_family_transformedGompertz:data_family_transformedWeibull                0.03
## fit_family_transformedNormal:data_family_transformedWeibull                 -0.20
## fit_family_transformedNormalIdentity:data_family_transformedWeibull         -0.14
## fit_family_transformedWeibull:data_family_transformedFréchet                -1.07
## fit_family_transformedBetaprime:data_family_transformedFréchet               0.49
## fit_family_transformedtrans.Normal:data_family_transformedFréchet           -0.16
## fit_family_transformedFréchet:data_family_transformedFréchet                 1.26
## fit_family_transformedGompertz:data_family_transformedFréchet               -0.86
## fit_family_transformedNormal:data_family_transformedFréchet                 -0.41
## fit_family_transformedNormalIdentity:data_family_transformedFréchet         -0.48
## fit_family_transformedWeibull:data_family_transformedBetaprime              -0.56
## fit_family_transformedBetaprime:data_family_transformedBetaprime             0.19
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime         -0.08
## fit_family_transformedFréchet:data_family_transformedBetaprime               0.25
## fit_family_transformedGompertz:data_family_transformedBetaprime             -0.58
## fit_family_transformedNormal:data_family_transformedBetaprime               -0.27
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime       -0.32
## fit_family_transformedWeibull:data_family_transformedGompertz                0.24
## fit_family_transformedBetaprime:data_family_transformedGompertz             -1.04
## fit_family_transformedtrans.Normal:data_family_transformedGompertz          -0.35
## fit_family_transformedFréchet:data_family_transformedGompertz               -0.94
## fit_family_transformedGompertz:data_family_transformedGompertz               0.48
## fit_family_transformedNormal:data_family_transformedGompertz                 0.06
## fit_family_transformedNormalIdentity:data_family_transformedGompertz         0.01
## fit_family_transformedWeibull:data_family_transformedtrans.Normal           -0.39
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal         -0.05
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal      -0.23
## fit_family_transformedFréchet:data_family_transformedtrans.Normal            0.17
## fit_family_transformedGompertz:data_family_transformedtrans.Normal          -0.27
## fit_family_transformedNormal:data_family_transformedtrans.Normal            -0.22
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal    -0.28
##                                                                          u-95% CI
## Intercept                                                                    2.08
## fit_link_transformedSoftplus                                                -0.13
## fit_family_transformedWeibull                                                0.24
## fit_family_transformedBetaprime                                              0.21
## fit_family_transformedtrans.Normal                                           0.50
## fit_family_transformedFréchet                                               -0.60
## fit_family_transformedGompertz                                              -0.07
## fit_family_transformedNormal                                                 0.18
## fit_family_transformedNormalIdentity                                         0.32
## data_link_transformedSoftplus                                                0.03
## data_family_transformedWeibull                                               0.23
## data_family_transformedFréchet                                               0.34
## data_family_transformedBetaprime                                             0.32
## data_family_transformedGompertz                                             -0.00
## data_family_transformedtrans.Normal                                         -0.24
## shapeHeavyTail                                                               0.05
## shapeRamp                                                                    0.10
## formulay~xPz1                                                                0.06
## formulay~xPz1Pz2Pz3                                                         -0.20
## fit_link_transformedSoftplus:data_link_transformedSoftplus                   0.25
## fit_family_transformedWeibull:data_family_transformedWeibull                 0.53
## fit_family_transformedBetaprime:data_family_transformedWeibull              -0.03
## fit_family_transformedtrans.Normal:data_family_transformedWeibull            0.26
## fit_family_transformedFréchet:data_family_transformedWeibull                -0.26
## fit_family_transformedGompertz:data_family_transformedWeibull                0.48
## fit_family_transformedNormal:data_family_transformedWeibull                  0.28
## fit_family_transformedNormalIdentity:data_family_transformedWeibull          0.34
## fit_family_transformedWeibull:data_family_transformedFréchet                -0.63
## fit_family_transformedBetaprime:data_family_transformedFréchet               1.03
## fit_family_transformedtrans.Normal:data_family_transformedFréchet            0.36
## fit_family_transformedFréchet:data_family_transformedFréchet                 1.76
## fit_family_transformedGompertz:data_family_transformedFréchet               -0.45
## fit_family_transformedNormal:data_family_transformedFréchet                  0.07
## fit_family_transformedNormalIdentity:data_family_transformedFréchet         -0.00
## fit_family_transformedWeibull:data_family_transformedBetaprime              -0.11
## fit_family_transformedBetaprime:data_family_transformedBetaprime             0.71
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime          0.46
## fit_family_transformedFréchet:data_family_transformedBetaprime               0.67
## fit_family_transformedGompertz:data_family_transformedBetaprime             -0.15
## fit_family_transformedNormal:data_family_transformedBetaprime                0.21
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime        0.16
## fit_family_transformedWeibull:data_family_transformedGompertz                0.72
## fit_family_transformedBetaprime:data_family_transformedGompertz             -0.61
## fit_family_transformedtrans.Normal:data_family_transformedGompertz           0.13
## fit_family_transformedFréchet:data_family_transformedGompertz               -0.57
## fit_family_transformedGompertz:data_family_transformedGompertz               0.95
## fit_family_transformedNormal:data_family_transformedGompertz                 0.54
## fit_family_transformedNormalIdentity:data_family_transformedGompertz         0.49
## fit_family_transformedWeibull:data_family_transformedtrans.Normal            0.04
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal          0.39
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal       0.24
## fit_family_transformedFréchet:data_family_transformedtrans.Normal            0.56
## fit_family_transformedGompertz:data_family_transformedtrans.Normal           0.14
## fit_family_transformedNormal:data_family_transformedtrans.Normal             0.22
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal     0.17
##                                                                          Rhat
## Intercept                                                                1.00
## fit_link_transformedSoftplus                                             1.00
## fit_family_transformedWeibull                                            1.00
## fit_family_transformedBetaprime                                          1.00
## fit_family_transformedtrans.Normal                                       1.00
## fit_family_transformedFréchet                                            1.00
## fit_family_transformedGompertz                                           1.00
## fit_family_transformedNormal                                             1.00
## fit_family_transformedNormalIdentity                                     1.00
## data_link_transformedSoftplus                                            1.00
## data_family_transformedWeibull                                           1.00
## data_family_transformedFréchet                                           1.00
## data_family_transformedBetaprime                                         1.00
## data_family_transformedGompertz                                          1.00
## data_family_transformedtrans.Normal                                      1.00
## shapeHeavyTail                                                           1.00
## shapeRamp                                                                1.00
## formulay~xPz1                                                            1.00
## formulay~xPz1Pz2Pz3                                                      1.00
## fit_link_transformedSoftplus:data_link_transformedSoftplus               1.00
## fit_family_transformedWeibull:data_family_transformedWeibull             1.00
## fit_family_transformedBetaprime:data_family_transformedWeibull           1.00
## fit_family_transformedtrans.Normal:data_family_transformedWeibull        1.00
## fit_family_transformedFréchet:data_family_transformedWeibull             1.00
## fit_family_transformedGompertz:data_family_transformedWeibull            1.00
## fit_family_transformedNormal:data_family_transformedWeibull              1.00
## fit_family_transformedNormalIdentity:data_family_transformedWeibull      1.00
## fit_family_transformedWeibull:data_family_transformedFréchet             1.00
## fit_family_transformedBetaprime:data_family_transformedFréchet           1.00
## fit_family_transformedtrans.Normal:data_family_transformedFréchet        1.00
## fit_family_transformedFréchet:data_family_transformedFréchet             1.00
## fit_family_transformedGompertz:data_family_transformedFréchet            1.00
## fit_family_transformedNormal:data_family_transformedFréchet              1.00
## fit_family_transformedNormalIdentity:data_family_transformedFréchet      1.00
## fit_family_transformedWeibull:data_family_transformedBetaprime           1.00
## fit_family_transformedBetaprime:data_family_transformedBetaprime         1.00
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime      1.00
## fit_family_transformedFréchet:data_family_transformedBetaprime           1.00
## fit_family_transformedGompertz:data_family_transformedBetaprime          1.00
## fit_family_transformedNormal:data_family_transformedBetaprime            1.00
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime    1.00
## fit_family_transformedWeibull:data_family_transformedGompertz            1.00
## fit_family_transformedBetaprime:data_family_transformedGompertz          1.00
## fit_family_transformedtrans.Normal:data_family_transformedGompertz       1.00
## fit_family_transformedFréchet:data_family_transformedGompertz            1.00
## fit_family_transformedGompertz:data_family_transformedGompertz           1.00
## fit_family_transformedNormal:data_family_transformedGompertz             1.00
## fit_family_transformedNormalIdentity:data_family_transformedGompertz     1.00
## fit_family_transformedWeibull:data_family_transformedtrans.Normal        1.00
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal      1.00
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal   1.00
## fit_family_transformedFréchet:data_family_transformedtrans.Normal        1.00
## fit_family_transformedGompertz:data_family_transformedtrans.Normal       1.00
## fit_family_transformedNormal:data_family_transformedtrans.Normal         1.00
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal 1.00
##                                                                          Bulk_ESS
## Intercept                                                                    7384
## fit_link_transformedSoftplus                                                21790
## fit_family_transformedWeibull                                                4657
## fit_family_transformedBetaprime                                              4906
## fit_family_transformedtrans.Normal                                           4441
## fit_family_transformedFréchet                                                4603
## fit_family_transformedGompertz                                               4574
## fit_family_transformedNormal                                                 4065
## fit_family_transformedNormalIdentity                                         3941
## data_link_transformedSoftplus                                               12494
## data_family_transformedWeibull                                               8466
## data_family_transformedFréchet                                               8607
## data_family_transformedBetaprime                                             8579
## data_family_transformedGompertz                                              8180
## data_family_transformedtrans.Normal                                          8195
## shapeHeavyTail                                                              14724
## shapeRamp                                                                   14513
## formulay~xPz1                                                               14542
## formulay~xPz1Pz2Pz3                                                         14188
## fit_link_transformedSoftplus:data_link_transformedSoftplus                  12532
## fit_family_transformedWeibull:data_family_transformedWeibull                 7124
## fit_family_transformedBetaprime:data_family_transformedWeibull               6807
## fit_family_transformedtrans.Normal:data_family_transformedWeibull            6346
## fit_family_transformedFréchet:data_family_transformedWeibull                 6061
## fit_family_transformedGompertz:data_family_transformedWeibull                6804
## fit_family_transformedNormal:data_family_transformedWeibull                  6322
## fit_family_transformedNormalIdentity:data_family_transformedWeibull          6180
## fit_family_transformedWeibull:data_family_transformedFréchet                 6133
## fit_family_transformedBetaprime:data_family_transformedFréchet               7679
## fit_family_transformedtrans.Normal:data_family_transformedFréchet            6557
## fit_family_transformedFréchet:data_family_transformedFréchet                 8506
## fit_family_transformedGompertz:data_family_transformedFréchet                6195
## fit_family_transformedNormal:data_family_transformedFréchet                  5794
## fit_family_transformedNormalIdentity:data_family_transformedFréchet          5763
## fit_family_transformedWeibull:data_family_transformedBetaprime               6816
## fit_family_transformedBetaprime:data_family_transformedBetaprime             7490
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime          6584
## fit_family_transformedFréchet:data_family_transformedBetaprime               6965
## fit_family_transformedGompertz:data_family_transformedBetaprime              6173
## fit_family_transformedNormal:data_family_transformedBetaprime                5845
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime        5800
## fit_family_transformedWeibull:data_family_transformedGompertz                6584
## fit_family_transformedBetaprime:data_family_transformedGompertz              6180
## fit_family_transformedtrans.Normal:data_family_transformedGompertz           5728
## fit_family_transformedFréchet:data_family_transformedGompertz                5659
## fit_family_transformedGompertz:data_family_transformedGompertz               6682
## fit_family_transformedNormal:data_family_transformedGompertz                 5793
## fit_family_transformedNormalIdentity:data_family_transformedGompertz         5753
## fit_family_transformedWeibull:data_family_transformedtrans.Normal            6194
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal          6704
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal       5629
## fit_family_transformedFréchet:data_family_transformedtrans.Normal            6547
## fit_family_transformedGompertz:data_family_transformedtrans.Normal           6221
## fit_family_transformedNormal:data_family_transformedtrans.Normal             5538
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal     5426
##                                                                          Tail_ESS
## Intercept                                                                    7596
## fit_link_transformedSoftplus                                                 6743
## fit_family_transformedWeibull                                                6266
## fit_family_transformedBetaprime                                              6607
## fit_family_transformedtrans.Normal                                           5942
## fit_family_transformedFréchet                                                6139
## fit_family_transformedGompertz                                               6723
## fit_family_transformedNormal                                                 5858
## fit_family_transformedNormalIdentity                                         5473
## data_link_transformedSoftplus                                                8009
## data_family_transformedWeibull                                               8234
## data_family_transformedFréchet                                               8317
## data_family_transformedBetaprime                                             8610
## data_family_transformedGompertz                                              7766
## data_family_transformedtrans.Normal                                          8376
## shapeHeavyTail                                                               8467
## shapeRamp                                                                    8341
## formulay~xPz1                                                                8046
## formulay~xPz1Pz2Pz3                                                          7669
## fit_link_transformedSoftplus:data_link_transformedSoftplus                   7717
## fit_family_transformedWeibull:data_family_transformedWeibull                 7559
## fit_family_transformedBetaprime:data_family_transformedWeibull               7822
## fit_family_transformedtrans.Normal:data_family_transformedWeibull            7340
## fit_family_transformedFréchet:data_family_transformedWeibull                 6523
## fit_family_transformedGompertz:data_family_transformedWeibull                7973
## fit_family_transformedNormal:data_family_transformedWeibull                  7334
## fit_family_transformedNormalIdentity:data_family_transformedWeibull          7683
## fit_family_transformedWeibull:data_family_transformedFréchet                 7432
## fit_family_transformedBetaprime:data_family_transformedFréchet               7719
## fit_family_transformedtrans.Normal:data_family_transformedFréchet            7617
## fit_family_transformedFréchet:data_family_transformedFréchet                 7418
## fit_family_transformedGompertz:data_family_transformedFréchet                7562
## fit_family_transformedNormal:data_family_transformedFréchet                  6664
## fit_family_transformedNormalIdentity:data_family_transformedFréchet          6805
## fit_family_transformedWeibull:data_family_transformedBetaprime               7940
## fit_family_transformedBetaprime:data_family_transformedBetaprime             7748
## fit_family_transformedtrans.Normal:data_family_transformedBetaprime          7344
## fit_family_transformedFréchet:data_family_transformedBetaprime               6998
## fit_family_transformedGompertz:data_family_transformedBetaprime              7533
## fit_family_transformedNormal:data_family_transformedBetaprime                7010
## fit_family_transformedNormalIdentity:data_family_transformedBetaprime        6974
## fit_family_transformedWeibull:data_family_transformedGompertz                7276
## fit_family_transformedBetaprime:data_family_transformedGompertz              7153
## fit_family_transformedtrans.Normal:data_family_transformedGompertz           7177
## fit_family_transformedFréchet:data_family_transformedGompertz                7340
## fit_family_transformedGompertz:data_family_transformedGompertz               7700
## fit_family_transformedNormal:data_family_transformedGompertz                 6951
## fit_family_transformedNormalIdentity:data_family_transformedGompertz         7054
## fit_family_transformedWeibull:data_family_transformedtrans.Normal            7304
## fit_family_transformedBetaprime:data_family_transformedtrans.Normal          8020
## fit_family_transformedtrans.Normal:data_family_transformedtrans.Normal       7629
## fit_family_transformedFréchet:data_family_transformedtrans.Normal            7461
## fit_family_transformedGompertz:data_family_transformedtrans.Normal           7605
## fit_family_transformedNormal:data_family_transformedtrans.Normal             6990
## fit_family_transformedNormalIdentity:data_family_transformedtrans.Normal     6783
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi    28.61      0.99    26.72    30.59 1.00    16296     8310
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

p1 = conditional_effects(m1, effects = "data_link_transformed:fit_link_transformed",
                    conditions = data.frame(
                      shape = NA,
                      data_family_transformed = NA,
                      fit_family_transformed = NA,
                      formula = NA)
                    )
plot(p1, plot = FALSE)[[1]] + scale_color_scico_d(palette = "batlow", end = 0.8) +
  xlab("Data Link") + ylab("AUC")+
  theme_bw(base_size = 12) +
  labs(color = "Fit Link", fill = "Fit Link")

ggsave(paste0(figure_path,"positive_auc_cond_link_link.pdf"), width = 210, height = (297/4)/1.5, units = "mm", useDingbats = TRUE)

p1$`data_link_transformed:fit_link_transformed` %>%
  ggplot(
    aes(y = estimate__,
        x = effect1__,
        ymin = lower__,
        ymax = upper__,
        color = effect2__)) +
  geom_pointrange(size = 0.4, linewidth = 0.8, position = position_dodge(width = 0.2)) +
  scale_color_scico_d(palette = "batlow", end = 0.8) +
  xlab("Data Link") +
  ylab("AUC ROC") +
  labs(color = "Fit Link") +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 9),
        strip.text.y = element_text(size = 9),
        axis.ticks.y = element_blank())

ggsave(paste0(figure_path,"positive_auc_cond_link_link.pdf"), width = 210, height = (297/4)*0.6, units = "mm", useDingbats = TRUE)

p2 = conditional_effects(m1, effects = "data_family_transformed:fit_family_transformed",
                    conditions = data.frame(
                      shape = NA,
                      data_link_transformed = NA,
                      fit_link_transformed = NA,
                      formula = NA))
plot(p2, plot = FALSE)[[1]] + scale_color_scico_d(palette = "batlow", end = 0.8) +
  xlab("Data Family") + ylab("AUC")+
   theme_bw(base_size = 12) +
        labs(color = "Fit Family", fill = "Fit Family")

ggsave(paste0(figure_path,"positive_auc_cond_fam_fam.pdf"), width = 210, height = (297/4)*1.2, units = "mm", useDingbats = TRUE)


p3 = conditional_effects(m1, effects = "fit_link_transformed",
                    conditions = data.frame(
                      shape = NA,
                      data_family_transformed = NA,
                      fit_family_transformed = NA,
                      data_link_transformed = NA,
                      formula = NA)
                    )
plot(p3, plot = FALSE)[[1]] + scale_color_scico_d(palette = "batlow", end = 0.8) +
  xlab("Fit Link") + ylab("AUC")+
   theme_bw(base_size = 12)

ggsave(paste0(figure_path,"positive_auc_cond_link.pdf"), width = 210/2, height = (297/4)*1.2, units = "mm", useDingbats = TRUE)


p4 = conditional_effects(m1, effects = "fit_family_transformed",
                    conditions = data.frame(
                      shape = NA,
                      data_family_transformed = NA,
                      fit_link_transformed = NA,
                      data_link_transformed = NA,
                      formula = NA)
                    )
plot(p4, plot = FALSE)[[1]] + scale_color_scico_d(palette = "batlow", end = 0.8) +
  xlab("Fit Family") + ylab("AUC")+
   theme_bw(base_size = 12)

ggsave(paste0(figure_path,"positive_auc_cond_fam.pdf"), width = 210, height = (297/4)*1.2, units = "mm", useDingbats = TRUE)


# Function to prepare a TPR/FPR DF averaging over different columns
roc_df <- function(cols){
  TPR <- calibration_df %>%
  filter(x_y_coef != 0) %>%
  group_by(across(all_of(cols))) %>%
  summarise(TPR_95 = mean(sig95),
            TPR_90 = mean(sig90),
            TPR_80 = mean(sig80),
            TPR_50 = mean(sig50),
            TPR_0 = 0,
            TPR_1 = 1)

FPR <- calibration_df %>%
  filter(x_y_coef == 0) %>%
  group_by(across(all_of(cols))) %>%
  summarise(FPR_95 = mean(sig95),
            FPR_90 = mean(sig90),
            FPR_80 = mean(sig80),
            FPR_50 = mean(sig50),
            FPR_0 = 0,
            FPR_1 = 1)

  df = merge(TPR, FPR)
  df = df %>%
    pivot_longer(
      cols = all_of(colnames(df)[! colnames(df) %in% cols]),
                    names_pattern = "TPR_(.*)$|FPR_(.*)$",
                    names_to = c("TPR", "FPR")) %>%
  mutate(TPR = as.numeric(TPR),
         FPR = as.numeric(FPR)) %>%
  mutate(CI = if_else(is.na(TPR), FPR, TPR),
         TPR = if_else(is.na(TPR), NA, value),
         FPR = if_else(is.na(FPR), NA, value))
  part1 <- filter(df, is.na(TPR))
  part1$TPR <- NULL
  part1$value <- NULL
  part2 <- filter(df, is.na(FPR))
  part2$FPR <- NULL
  part2$value <- NULL
  df = merge(part1, part2)
  
  return(df)
}

# Full ROC Plot
df <- roc_df(cols = c("formula", "fit_link_transformed", "fit_family_transformed", "shape", "data_link_transformed", "data_family_transformed"))

## `summarise()` has grouped output by 'formula', 'fit_link_transformed',
## 'fit_family_transformed', 'shape', 'data_link_transformed'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'formula', 'fit_link_transformed',
## 'fit_family_transformed', 'shape', 'data_link_transformed'. You can override
## using the `.groups` argument.

roc <- df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1),labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_family_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_family_transformed), data = filter(df, CI == "95"), size = 2) +
  facet_grid(shape + data_link_transformed + data_family_transformed ~ formula + fit_link_transformed) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc

ggsave(paste0(figure_path, "positive_roc.pdf"), width = 210*2, height = (297/4)*16, units = "mm", useDingbats = TRUE)

# Compare the fit configurations and average over all DGPS
df <- roc_df(cols = c("formula", "fit_link_transformed", "fit_family_transformed"))

## `summarise()` has grouped output by 'formula', 'fit_link_transformed'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'formula', 'fit_link_transformed'. You can
## override using the `.groups` argument.

roc_by_fit <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_family_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_family_transformed), data = filter(df, CI == "95"), size = 2) +
  facet_grid(fit_link_transformed ~ formula) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_fit

ggsave(paste0(figure_path, "positive_roc_by_fit.pdf"), width = 210, height = (297/4) * 1.5, units = "mm", useDingbats = TRUE)

# Compare the DGPs and average over all fit configurations
df <- roc_df(cols = c("shape", "data_link_transformed", "data_family_transformed", "fit_family_transformed"))

## `summarise()` has grouped output by 'shape', 'data_link_transformed',
## 'data_family_transformed'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'shape', 'data_link_transformed',
## 'data_family_transformed'. You can override using the `.groups` argument.

roc_by_dgp <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_family_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_family_transformed), data = filter(df, CI == "95"), size = 2) +
  facet_grid(shape + data_link_transformed ~ data_family_transformed) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_dgp

ggsave(paste0(figure_path, "positive_roc_by_dgp.pdf"), width = 210, height = (297/4) * 3, units = "mm", useDingbats = TRUE)

# Compare the DGP shape and fit configurations 
df <- roc_df(cols = c("fit_link_transformed", "fit_family_transformed", "shape"))

## `summarise()` has grouped output by 'fit_link_transformed',
## 'fit_family_transformed'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'fit_link_transformed',
## 'fit_family_transformed'. You can override using the `.groups` argument.

roc_by_shape_and_fit <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_family_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_family_transformed),
             data = filter(df, CI == "95"), size = 2) +
  #geom_point(aes(color = fit_family_transformed), alpha = 0.5) +
  facet_grid(fit_link_transformed ~ shape) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_shape_and_fit

ggsave(paste0(figure_path, "positive_roc_shape_and_fit.pdf"), width = 210, height = (297/4) * 1.2, units = "mm", useDingbats = TRUE)

# Compare the DGP Likelihood and fit configurations 
df <- roc_df(cols = c("fit_link_transformed", "fit_family_transformed", "data_family_transformed"))

## `summarise()` has grouped output by 'fit_link_transformed',
## 'fit_family_transformed'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'fit_link_transformed',
## 'fit_family_transformed'. You can override using the `.groups` argument.

roc_by_shape_and_fit <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_family_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_family_transformed),
             data = filter(df, CI == "95"), size = 2) +
  #geom_point(aes(color = fit_family_transformed), alpha = 0.5) +
  facet_grid(fit_link_transformed ~ data_family_transformed) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_shape_and_fit

ggsave(paste0(figure_path, "positive_roc_dgp_lik_and_fit.pdf"), width = 210, height = (297/4) * 1.2, units = "mm", useDingbats = TRUE)

# Compare the links
df <- roc_df(cols = c("fit_link_transformed", "data_link_transformed"))

## `summarise()` has grouped output by 'fit_link_transformed'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'fit_link_transformed'. You can override
## using the `.groups` argument.

roc_by_shape_and_fit <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_link_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_link_transformed),
             data = filter(df, CI == "95"), size = 2) +
  #geom_point(aes(color = fit_family_transformed), alpha = 0.5) +
  facet_grid(. ~ data_link_transformed) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_shape_and_fit

ggsave(paste0(figure_path, "positive_roc_links.png"), width = 210*2, height = (297/4) * 1.2*2, units = "mm")#, useDingbats = TRUE)

# Compare the likelihoods
df <- roc_df(cols = c("fit_family_transformed", "data_family_transformed"))

## `summarise()` has grouped output by 'fit_family_transformed'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'fit_family_transformed'. You can override
## using the `.groups` argument.

roc_by_shape_and_fit <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_family_transformed), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_family_transformed),
             data = filter(df, CI == "95"), size = 2) +
  #geom_point(aes(color = fit_family_transformed), alpha = 0.5) +
  facet_grid(. ~ data_family_transformed) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=c(0, 1, 8, 4, 5, 17, 6, 7)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_shape_and_fit

ggsave(paste0(figure_path, "positive_roc_likelihoods.png"), width = 210*2, height = (297/4) * 1.2*2, units = "mm")#, useDingbats = TRUE)

# Compare LLs
df <- roc_df(cols = c("fit_family_transformed", "fit_link_transformed", "data_family_transformed", "data_link_transformed")) %>%
  mutate(data_LL = paste0(data_family_transformed, data_link_transformed),
         fit_LL = paste0(fit_family_transformed, fit_link_transformed))

## `summarise()` has grouped output by 'fit_family_transformed',
## 'fit_link_transformed', 'data_family_transformed'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'fit_family_transformed',
## 'fit_link_transformed', 'data_family_transformed'. You can override using the
## `.groups` argument.

roc_by_LL <-  df %>%
  ggplot(aes(FPR, TPR)) +
  coord_cartesian(xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_continuous(breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  geom_line(aes(color = fit_LL), linetype = "longdash", linewidth = 1) +
  geom_abline(color = "grey", linetype = "dotted", slope = 1) +
  geom_point(aes(FPR, TPR, shape = fit_LL),
             data = filter(df, CI == "95"), size = 2) +
  #geom_point(aes(color = fit_family_transformed), alpha = 0.5) +
  facet_wrap(. ~ data_LL) +
  scale_colour_scico_d(palette = "batlow", end = 0.8) +
  scale_shape_manual(values=seq(1,16)) +
  theme_bw(base_size = 12) +
  theme(strip.text.x = element_text(size = 10),
        strip.text.y = element_text(size = 10),
        legend.position = "bottom") +
  labs(color = "Fit family", shape = "Fit family") +
  xlab("False positive rate") + ylab("True positive rate")

roc_by_LL

ggsave(paste0(figure_path, "positive_roc_LL.png"), width = 210*2, height = (297/4) * 1.2*3, units = "mm")#, useDingbats = TRUE)